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In quasi-static nonlinear time-dependent analysis, the choice of the time discretization is a complex 
issue. The most basic strategy consists in determining a value of the load increment that ensures 
the convergence of the solution with respect to time on the base of preliminary simulations. In more 
advanced applications, the load increments can be controlled for instance by prescribing the number of 
iterations of the nonlinear resolution procedure, or by using an arc-length algorithm. These techniques 
usually introduce a parameter whose correct value is not easy to obtain. In this paper, an alternative 
procedure is proposed. It is based on the continuous control of the residual of the reference problem 
over time, whose measure is easy to interpret. This idea is applied in the framework of a multiscale 
domain decomposition strategy in order to perform 3D delamination analysis. 

1 INTRODUCTION 

The virtual testing of delamination is an objective widely spread among industrialists especially in the 
aeronautical field. To achieve it, two research thematics which have undergone large evolution during 
the last twenty years need to be put in conjunction: the pertinent modeling of composites and the 
efficient computation of structures. 

Indeed, there have been many advances toward a better understanding of the mechanics of lam- 
inated composites and of damage mechanisms. Two kinds of modeling have proved their validity: 
microscale and mesoscale models. Microscale models are strongly connected to the physics of the ma- 
terial and thus provide a reliable framework for simulation |12[ I22j . Unfortunately, the computation 
of models defined at the micro scale require such a fine discretization that only small test specimens 
can be simulated, structural computations being out of reach even on recent hardware. Meso-models 

m [13 [7j are defined at a scale which enables both the introduction of physics-based ingredients 
and the simulation of small industrial structures. They indeed most often rely on the definition of two 
meso-constituents, the ply (3D entity) and the interface (2D entity), which are modeled using contin- 
uum (damage) mechanics, their behavior being obtainable from the homogenization of micro-models 
|21j . Anyhow, for reliable simulation, discretizations still need to be fine (in order, for instance, to 
represent correctly the gradients of stresses due to edge effects which are responsible for the initiation 
of many degradations) and associated systems thus remain very large (in terms of number of degrees 
of freedom) and strongly nonlinear (with potential instabilities) . 

As a first approach of the reliable simulation of quasi-static simulations of the delamination in 
composite structures, we chose in [TU] to neglect the effect of deterioration within the plies and to 
lump the degradations in the interfaces. We thus retained the mesomodel presented in |5] where 
the delamination ability is localized in the interfaces and handled through a cohesive behavior. The 
space discretization is considered sufficiently fine to represent accurately any evolution of multiple 
delamination cracks (sufficient number of Gauss points in the length of the process zone [27l[Tl[7]). At 
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each time step of an incremental time discretization scheme, the associated large nonlinear system is 
solved using a three-scale domain decomposition strategy. Based on the mixed LaTIn-based domain 
decomposition method |14| , this strategy has been given high numerical efficiency by adapting various 
ideas from the work of [23l|24l[25] to the computation of delamination. Three-dimensional simulations 
of the delamination in realistic structural components have been performed on parallel computers 
without the need to perform local space refinement. 

Though, a complex issue arises when choosing the load increments: the solution to softening quasi- 
static problems depends on the time discretization scheme parameters (non-uniqueness of the solution 
and possible bifurcation paths). This remark brings us to the field of the validation. In the literature, 
numerous error indicators have been developed to control a posteriori the global error introduced 
in finite element schemes for linear problems [51 [551 [20] • These indicators have been extended to the 
validation of nonlinear time-dependent problems [T31[S1[S|. One of the most advanced criterion is the so- 
called error in the constitutive law [13| . A solution to the nonlinear evolution problem being computed 
using a FE scheme and a classical time integration procedure, one constructs a solution which satisfies 
the kinematic and static admissibilities, and lump the residual of the nonlinear evolution problem 
equations in the constitutive laws. A measure of this residual permits to control at the same time 
the discretization error in space, in time and the error introduced by the iterative solution procedures 
[IB] . This idea has been formalized in [TT] for materials described using internal variables. The 
state equations are satisfied by the reconstructed solutions, the measure of the non-verification of 
the evolution laws permits to derive a strict upper bound to the solution error. Though, this new 
admissible solution is not easily constructed in the case of softening behaviors. Specific developments 
in |17| meant to tackle this difficulty, and the resulting procedure is used in jl6j to derive an adaptive 
refinement procedure in space and time. Note that, at the present time, a link between the error in the 
constitutive law and the error in the solution is still to be established in the case of softening materials. 

The aim of the work presented in this paper is double. The first is to define a comprehensive time 
discretization error indicator inspired from the work of |13( I16j for delamination analysis and to ensure 
that its computation and use is numerically efficient within the LaTIn-based domain decomposition 
strategy. Our second goal is to use the developed indicator to control on the fly the load increment in 
quasi-static analysis in order to ensure the convergence of the computed solution. 

The paper is organized as follows. The reference delamination problem is presented in Section [2| 
The dependency of the solution to this problem on the time discretization scheme is demonstrated. In 
the following section, we present a time-dependent error indicator based on the error in the constitutive 
law and computed with respect to a continuous solution in time, constructed by interpolation over 
each time step. Although very general, this indicator is not directly suitable for the LaTIn-based 
multiscale strategy used to perform the nonlinear resolutions. The main features of this strategy are 
recalled in Section [4j We focus in particular on the indicator based on the error in the constitutive law 
used to estimate the convergence of the iterative procedure. In Section [sj this convergence indicator 
is associated to the previous developments to derive an alternative and cheap time discretization error 
indicator, which is the basis for the development of an automatic time-step-control procedure. At 
last, this technique is validated on multiscale and parallel delamination simulations in Section [6j Two 
different problems are assessed: a simple and stable problem in which the time increments correspond 
to the increases in the prescribed load, and a more complex and unstable problem, solved using an 
arc-length procedure, in which the time increments correspond to the value of the arc-lengths. 

2 THE REFERENCE PROBLEM AND ITS DISCRETIZA- 
TION IN TIME 

2.1 Reference problem at a given time of the analysis 

The delamination simulation is performed under the assumptions of quasi-static, isothermal evolution 
over time and small perturbations. 
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The laminate structure E occupying Domain f2 is made out of Np adjacent plies occupying 
Domains (i^p)pg|i, Nr} (^^ boundaries {dflp)p^^i_ Np})j separated by (Np — 1) cohesive interfaces 
(^p)peli. Np-ij and (see Figure ([6|, Page 10 1. An external traction field (respectively a displace- 



ment field U_^) is applied to the structure on Part dilf (respectively dflu) of the boundary dQ of 
Domain ft. The volume force is denoted /^. Let Up be the displacement field, the Cauchy stress 
tensor and the symmetric part of the displacement gradient in Ply P. 

At every time t € [0 T] of the analysis, the reference non-linear equilibrium problem reads: 

Find Sref — (sp)pe|i, Np}), where sp = {up,g^p), which satisfies the following equations: 
• Kinematic admissibility on dflu'- 



Global equilibrium of Structure E: V(up*)pg|j^ 

J2 f Tv L eSup*)) dn f f^.Up*dQ-J2f F^.Up*dT 

p Jup ^ ^ p Jnp p Janpnafif 



(1) 



(2) 



where [u]^ is the jump of displacement of Interface Ip: [u]^ = Up^i — Up and np is the outer 
normal to the boundary dQ,p. 

• Linear orthotropic behavior of the plies: 



(3) 



Constitutive law of the cohesive interfaces, local on any interface Ip. The elastic damageable 
law proposed in 0] is described using continuum damage mechanics. Three internal variables 
(rfi)ie|i, 3] (one for each delamination mode: traction along np and shear along ti and ^2 oii 
Figure (U])), ranging from to 1 are introduced in the surface strain energy in order to take 
into account the irreversible damage mechanisms. 






np 










Figure 1: The mesomodel entities 



Two state equations are derived from the expression of the free energy. The first one 
establishes a relation between the dual interface unknown g^pTip, and the primal interface 
unknown \u\ „ : 



gp-np 



— which gives ^p.np=Ip^MpJ|,,[o jMp 
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where, in the basis (rip, fj^, ^2)1 ^+ being the positive indicator function: 



|re[0 i] 



(l-di)fc.° 











(1 - d2)k° 



The second state equation Unks the thermodynamic forces (Yi)i^^i^ 3| to the primal interface 
unknown: 

1 „ / \2 



Y,, 



ddi 



where 



Y2 
Y^ 



-fcO (h+{[v\^.np) 



(5) 



The evolution laws are: 

di — d2 = ^ min{l,w(y)} 



where < 



(6) 



Further details on this cohesive zone model and identification issues can be found in |4j. 

The dissipated energy Edissi will be used in this paper as a global measure of the delaminated 
area of the cohesive interfaces: 



Edi 



dissi 



p Jip Jo \^^^ j p Jip 



AddT 



(7) 



where A is a scalar which only depends on the parameters of the interface model. 



In the following developments, the investigations are restricted to simulations under prescribed 
forces and displacements following a unique load function of time. In this context, the volume force 
will be assumed negligible. These assumptions are not mandatory to make use of the work presented 
in this paper, but simplify the construction of a continuous solution over time (Section [s]). 



2.2 Time discretization scheme 

An incremental procedure is used to solve the problem over time. It consists in discretizing the time 
of the analysis [0 T] in A'^ intervals [tn t„+i]„g|o Ar-i|. Successive nonlinear problems are solved at 
each computation time (tn)rae[o, n}- 

Hence, a solution to the discretized problem in time is a set of A^ + 1 solutions satisfying the 
reference problem equations, the time dependency in the constitutive laws being discretized. More 
precisely, at Computation time t„+i, the discretization of Equations Q and ^ reads: 




In general, the time discretization is chosen so that within each interval [t„ t„-|_i]„g|o Ar_i|, the 
evolution of the prescribed load is monotonic, which will also be assumed in the following. 
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2.3 Influence of the time increments on the solution to the discretized 
delamination problem 

The solution to the discretized reference problem reached at time T strongly depends on the time 
increments for two reasons: 

• the discretized cohesive law (Equation ([s])) depends on the discrete history of the interface 
variables. Hence, the residual stiffness of the cohesive interfaces depends on the time increments. 
This phenomenon is illustrated in the next section. 

• structural problems involving softening materials may be unstable and may have multiple solu- 
tions. In those cases, the solution paths depend on both the algorithm used at each computation 
time step and the initialization of this algorithm {i.e.: the previous converged solution). The 
resulting dependency on the time increments will be demonstrated in the last section of this 
paper. 



DCB (double cantilever beam) test case The laminate structure that we consider is made out 
of four isotropic plies (Figure [s]) . One part of the median cohesive interface is replaced by a contact 
interface in order to simulate an initial crack in the structure. Displacements are prescribed for the 
crack to propagate in a stable manner. The final prescribed displacement is set to a predefined value, 
which fixes the propagation length. The initial stiffness of the cohesive interfaces is obtained by 
integrating the Young and shear moduli of the matrix in the "thickness" of the interface (1/10 the 
thickness of the plies) [3] . 

The solution is not unique and depends on the load increments. Figure [3] presents the damage state 
in the upper cohesive interface, four different time discretizations being applied (these results will be 
fully detailed later on, for the values of the successive load increments are obtained by the adaptive 
time step procedure described in Section [5|. j^*™^''''' jg the criterion driving the time discretization 
(the largest j/*™'^''^''^ the coarser the discretization). In cases 1 and 2, the number of time increment 
used in the propagation phase of the analysis are, respectively, 69 and 21. The differences in the 
damage state of the interfaces are not significant, the evolution of the crack front being sufficiently 
slow to capture the effects of the stress concentrations. Hence, both these solutions are converged 
with respect to the time. In case 3, obtained with 9 coarse time increments, the solution is slightly 
different from the previous reference cases. Finally, in case 4, using only 5 time increments to describe 
the propagation of the crack clearly leads to the appearance of damage strips in the upper and lower 
interfaces. This is due to the effect of the stress concentration at the tip of the crack which propagates 
in a discrete manner with respect to time. 




Figure 2: Definition of the four-ply DCB problem 
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Figure 3: Influence of the prescribed value jy**™'^''''^ on the damage state in the upper cohesive interface 
of the DCB problem 

3 A TIME DISCRETIZATION ERROR INDICATOR 

We suppose that two consecutive solutions to the reference problem, s„ at Time t„ and s„+i at Time 
have been computed using a nonlinear resolution strategy. The aim is to evaluate the relevancy 
of the solution computed at Time t„+i, the continuous evolution of the structure over the current time 
step [tn tn+i] being a priori unknown. We propose to construct an interpolated solution over the time 
step in order to monitor the residual of the nonlinear reference problem continuously. 

3.1 Interpolation of the kinematic and static fields over a time step 

Let us prescribe the continuous evolution of the prescribed boundary values over the time step: 



Vt € [tn tn+l], 



VAf e aa., C/,|, = a(t) C/,|,^ + (1 - a{t))U,^^^^^ 



(9) 



where the function a{t) is the restriction of the load function over [i„ in+i]- In the case of a linear 
evolution (which will be the case in our applications), it simply reads: 



\ft e [tn tn+l], a{t) 



t tri 



(10) 



The evolution of the kinematic and static fields over the current time is assumed to follow the 
evolution of the prescribed loading (see Figure Q), which writes: 



V<e [tn VPe II, Npj, 



up\t = "(^Mp|t„ + (1 - a(t))up|t^^^ 



(11) 



Sn and Sn+i are two solutions of the reference problem. In particular, they satisfy the following set 
of linear equations: 

• kinematic admissibility. Equation ([I]) 
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Figure 4: Schematic representation of the interpolation performed over each time step 

• static admissibihty, Equation ([2|, the volume force being assumed negligible. 

• constitutive law of the plies, Equation ^ 

As a consequence, the interpolated kinematic and static fields over the current time step also satisfy 
this set of linear equations. Hence, the residual of the reference problem at any time t € [tn tn+i] 
is the residual of the constitutive law of the cohesive interfaces, which remains the only non-satisfied 
equation. 



3.2 Evolution of the damage variables over the current time step 

At any intermediate time t e [tn in+i], the internal variables are calculated with respect to the 
continuous history of the interpolated displacement field on Time interval [0 t\ . Let us define a new 
stress field d_ which satisfies the nonlinear constitutive law of the interfaces: 

Vie [t„ t„+i],VPe II, iVp-il, on/p. 




Alternatively, one can update the damage variables with respect to the interpolated stress field, 
and define a jump of displacement field [u] satisfying the constitutive law of the cohesive interfaces. 

The damage variables initially computed at time t„+i by the nonlinear resolution strategy are 
discarded. Indeed, they may differ from the ones obtained at time tn+i by the continuous construction 
over [tn tn+i], for solution s„+i only satisfies the discretized cohesive law ([8|. The residual of the 
reference problem equations at Time tn+i obtained when updating the damage variables can be reduced 
by lowering the time increment At ~ tn+i — tn and performing new nonlinear resolutions at Time tn+i, 
which will be detailed in Section O 



3.3 Definition of the time discretization error indicator 

A measure ["interp" stands for "interpolation") of the residual of the reference problem equa- 

tions at any time i g [i„ tn+i] can be obtained by summing the local contributions to the error in the 
nonlinear constitutive laws: 

z^lf-^^^ " ^^^^^ I, ^ 1,^^^ /■ ^T^^p (13) 

P II (^p|f + |p|f)«p ll/P 
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Or alternatively if the history is updated with respect to the interpolated stress field 



v'^^'^^'p = J2 f^^' zz^^'i (14) 



u\ + u, 




Figure 5: Schematic representation of the time discretization error indicator 



The time discretization error indicator at Time tn+i is defined as the maximum value of the previous 
measure over [t„ (see Figure (Is])), which reads: 



i^fr' = max or ahernatively Jf/'"^ = max J^,"'*""^ (15) 

The concept introduced here finds its roots in the work of [131 [5] , in which the sum over time of the 



product of Criteria ( 13 1 and ( 14 ) is used to measure the error in the constitutive law due to both space 
and time discretization for materials satisfying Drucker's stability equality. Three main differences 
should be outlined here: 

• In the case of softening materials, Drucker's stability equality is not satisfied. The mathematical 
properties which result from the definition of the Drucker's law-based criterion do not apply. 
Hence, making use of this criterion is not relevant. In addition, computing j/*™'^ requires the 
monotony of the interface behavior (uniqueness of the admissible displacement jump for any 
arbitrary local stress state). In the following developments, we will use Criterion ^y*""'^ to measure 
the residual of the reference problem equations over the current time step. 

• Our final goal being to provide an algorithm to control "on-the-fly" the time increments, i/*™^ 
is not a norm over the whole time of the analysis, but it instead is evaluated locally over each 
time increment. 

• To be consistent with [T51I5] the field should also be reconstructed with respect to the space 

variables so that it satisfies exactly the static admissibility condition (pi). In this work we focus 
on the time discretization and so we content ourselves with a weak (discrete) static admissibility. 
At Times t„ and tn+i solution fields satisfy the constitutive law of the plies ([s]), the kinematic 
admissibility and the static admissibility "in the finite element sense". Thus Criterion ^y**™^ 
(which is introduced without reference to space discretization) only accounts for the error due 
to time discretization. 
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3.4 Practical considerations 



3.4.1 Sub-intervals 



In practice, jy*"*^'"? is computed at a given set of intermediate times within the current time step. 
[tn tn+i] is subdivided into Ng subintervals [ti i^+iJjgiQ^ Ws-il: the time discretization error indicator 

,iir 



lyttme |-,gijjg computed as: 



1/ ™" = max (16) 



3.4.2 Error in the cohesive law 

Computing requires to extract the transverse constraints Np} which is not directly 

available in finite element codes. Usually, cohesive interface elements are used to overcome this problem. 
Classical incremental Newton solvers can then be used to solve the delamination problem at each 
computation time (in)ne|o, nj- The technique to control the time increment that we propose in 
Section [5.1.2| can directly be applied to such approaches. 

We focus on the insertion of the control technique within the framework described in [TD]. The 
principle is to use an incremental LaTIn-based domain decomposition strategy |18j to efficiently solve 
(in parallel) the delamination problem at each computation time. In this case, the cohesive behavior is 
directly described as a nonlinear joint between substructures. The mixed description of the interface 
behavior makes the transverse constraints available naturally. As it shall be detailed in Section [5] 
the time discretization error indicator can be defined as a time-dependent version of the convergence 
indicator used to stop the iterations of the LaTIn algorithm. 



4 THE NONLINEAR RESOLUTION STRATEGY 

We propose an overview of the domain decomposition strategy used to perform the successive nonlinear 
resolutions of the delamination analysis, first in the stable case, then in the unstable case, where it is 
combined with an arc-length procedure. We focus in a second time on the development of a convergence 
indicator based on the error in the constitutive law |13| to stop both of these iterative solvers. Further 
details concerning the multiscale and parallel computing aspects can be found in jlOj . 



4.1 Substructured formulation of the reference problem 

The laminate structure E is decomposed into substructures and interfaces as represented in Figure 
^ . Each of these mechanical entities possesses its own kinematic and static unknown fields linked by 
its behavior. The substructuring is driven by the will to match domain decomposition interfaces with 
material cohesive interfaces, so that each substructure belongs to a unique ply and has a constant linear 
behavior. Each substructure is defined in a domain 51^ such that E G |1, ueJ {jie being the total 
number of substructures) and is connected to a adjacent substructures through interfaces Tee' = 
dVlE r\ dflE' where E' e |1, n^J (Figure Q). The surface entity Tee' applies force distributions 
F^^, F^, as well as displacement distributions W^, W^, to Substructure E and Substructure E' 
respectively. On Substructure E such that c^JIb n dfl ^ 0, the boundary condition {Uj^,Fj) is applied 
through a boundary interface r^;^ . Let us deHne T e = [J e' eli nsl ^sb'UF^;^. We finally introduce a^, 
the Cauchy stress tensor, and e(u^), the symmetric part of the displacement gradient, in substructure 
E. 

The substructured quasi-static problem at any computation time of the time discretization 
scheme consists in finding s = {sE)Eeli, where se = (W F^), which satisfies the following 
equations: 

• Kinematic admissibility of Substructure E: 

Uew. = We (17) 
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Figure 6: Substructuring of the laminated composite structure 

• Static admissibility of Substructure E: V(u^^, W^*) eUs ^ / 'UE*\dnE ~ ^e*^ 

• Linear orthotropic behavior of Substructure E: 

• Behavior of the interfaces Tee' S ^e'- 

nEE'{WE.WE.,FE,FE,)^Q 

• Behavior of the interfaces F^;^ G (F^ n dVl): 

■ReA^e^Fe) = {We = Md on dflu and = F^ on d^f) 

In delamination analysis, the formal relation TZee' — ^ reads: 





For perfect interface: 



• For cohesive interface: 



Fe + Fe, 



W, 



We, = 



Fe + Fe, 







FE=Kp{iWE,-WE)^,^^,„^ 



{We, -We) 



(18) 



(19) 



(20) 



(21) 



where Substructure E (respectively E') belongs to Ply P (respectively P +1). 

4.2 Iterative resolution of the stable nonlinear substructured problem 

The equations of the problem can be split into the set of linear equations in substructures (static 
and kinematic admissibility of the substructures, linear constitutive law of the substructures) and 
the set of local equations in interface variables (behavior of the interfaces). The solutions s = 
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(s_e)£;g[i, tie] = (Vt^g. F^) Kfzp „g| to the first set of equations belong to Space Ad, while the so- 
lutions s = (s£;)_Bg|i, „e1 — ^.F^)fp,(zp „g| to the second set of equations belong to T. Hence, the 
converged solution Sref is such that Sref € Ad f] T. 

The LaTIn resolution scheme consists in searching for the solution s^e/ alternatively in these two 
spaces along search directions E+ and E~ (see Fig. [8|: 

• Find s*+3 g r such that ^s*+3 — s*^ G E+ (local stage) 

• Find s'+^ e Ad such that ^s'+^ — s*+3^ e (linear stage) 
In the following, the subscript i will be dropped. 

Local stage One searches for a solution s = (F ^) p,(zp satisfying the local equations on 
the interfaces {TZee' = or TZe^ =0), and search direction equation E+, introduced locally on the 
interfaces : 

(Ie ~ Fe) - k+i^E ~We) = G (22) 

At this stage, variables F^e W-e known from the previous semi-iteration. 

In the case where TZee' = is a nonlinear equation, the local problem is solved by a quasi-Newton 
algorithm. 
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Linear stage One searches for a solution s = (F^, W^^) p,^p verifying the linear equations on 
each substructure and, at best, a search direction equation local on the interfaces, under the 
constraint of average equihbrium of the interface forces : 

I F^ip^ = argmin |^ (^(^s " ^E? + {Ee " Ie)-^^ - g^)) 
I under the constraint: y{E',E), IlM Fr, +nM Fr,, =0 

The macroscopic projectors Ilj'ji^ ^ extract an average of the interface forces, which is transfered into 
the whole structure. Technicahy, this stage consists in solving, in parallel, independent linear problems 
on the sub-structures (using finite elements) and a small macroscopic linear problem which is global 
over the structure (and discrete by construction). 



(23) 



4.3 Iterative resolution of the unstable nonlinear problem 

When a snap-back appears in the global behavior of the simulated structure, the incremental LaTin 
algorithm is switched to a well-known local arc-length algorithm [37| [51 [TU] . The algebraic nonlinear 
problem to solve at Time in an unstable phase, reads: 

9int (f^|t„ + i, (t^|T)r<t„+i) - A|t„^^ goxt = (24) 

The amplitude of the loading \\t„j_i is unknown. A control equation inspired from [2 is introduced so 
that the maximum local increment in the jump of displacement over all the cohesive interfaces takes 
a predefined value A^: 

c(Af/|,„^JAC/|,„^, =AZ (25) 

where the A . unknowns are the increments in the quantities over Time step [t„ in+i]- c is then the 
operator which extract the maximum of the (positive) jump increment. 

Classically, the non-linear system (24 25 ) is solved by a modified Newton- Raphson scheme: 



• The linearization of (24) and (25) around point (C/',A*) leads to the system to solve at the 
prediction step of the (i -|- l)"* iteration of this scheme: 

A/ + c(A^|-^^j£/|,„ 



c(AC/|^^_^^Jk(;7|',^^^,(C/|,),<,„^,) \ext (26) 



The inversion of the linearized stiffness operator (i.e.: the resolution of the linear system U — 
K(i7|'j (?7|r)T<t„+i )~^9cxt) is performed by using the domain decomposition method described 
previously (the internal variables of the interfaces are fixed during the resolution) 

• The correction step of the Newton scheme consists in updating Operators K and c with respect 
to the kinematic field U^^^^^ found at the prediction stage. 



4.4 Stopping criterion 

4.4.1 Stable phase (LaTIn algorithm) 

In order to evaluate the convergence of the LaTIn algorithm, one classically measures the distance 
between spaces Ad and T along search direction E~ [11] (criterion labeled i^"'' on Figure "sd" 
standing for "search direction" ) . In the work of |10[ I19j , a new criterion based on the error in the 
constitutive law has been successfully assessed (in order, originally, to get rid of the dependency of 
Convergence indicator t^^"* on the parameters of the LaTIn solver). The solutions resulting from a 
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Figure 9: Classical convergence indicator i/"'^ of the LaTIn solver and indicator v'^*'^^ based on the error 
in the constitutive law 



linear stage of the LaTIn algorithm satisfy all the equations of the substructured reference problem 
except the interfaces laws (201 and (21 1, whose residuals can be easily computed (Figure More 
precisely, a solution s*^^ € Ad being reached, an indicator of the convergence of the algorithm is given 
by integrating the corresponding local residuals of the interface behaviors over the structure (residuals 
of Equations ^ and ^ evaluated for s'+^ = is'^^)E<^li, n^} = iW!+\ F!+^)Eeii, n^})- 

Let us call q the number of interface relations being used (i.e.: the number of distinct interface 
behaviors (JZee' — ^)[E.E')£li, ^nd [TZe^ — ^)e£Ii, ueI used in the structure). In our case, q = A 
(perfect interfaces, cohesive interfaces with homogeneous constants, Dirichlet and Neumann boundary 
conditions). Ti is the set interfaces of Behavior i, for all i S |1, q\. The vectorial relations TZee' — 
for i £ |1, q} and Tee' & or F^^ € f ■ are made of pi vectorial equations Qij = (2 equations for 
cohesive or perfect interfaces in 3D, 1 equation for boundary interfaces). Here, subscript j ranges from 
1 top. Convergence indicator {^'Hter" stands for "iterative") reads: 
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where 



,2 reTi 



i=l 3 = 1 



pgr, ^ 



(27) 



where, in the case of delamination (i.e : involving perfect and cohesive LaTIn interfaces): 
• on a perfect interface F^;^;' G Fi : 



Qii = - Fe, 



Qi2 = We- We, 

Qi2^We + We, 



(28) 



on a cohesive interface Tee' € r2 '■ 



Q21 =Fe-Kp {{We, - WE)ite{t„^„[o *„]}) {We, ~ We) 
Q21 =Fe+Kp {{We, - WE)ite{t„^,.^[o t„]}) {We, - We) 
Q22 = Fe, - Kp {{We, - WE)it^{t„^,,[o *„]}) {We - We,) 
Q22 = Fe, + Kp {{We, ~ lEB)|t6{t„^„[o *„]}) {We - We,) 



(29) 



where P G |1, Np — 1]. Note that, in Equation (29 1, the history of the interface variables 



during the current load increment is not taken into account, for it is unknown at this stage of 
the resolution procedure. 
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• on an interface transmitting Neumann's boundary condition : 

Ssi-Zij-Zd Q3i=ZB + Zd (30) 

• on an interface transmitting Dirichlet's boundary condition Te^ : 

Q41 ^We-W^ Q41 ^We + W^ (31) 

The computation of this criterion is very cheap as it simply requires local integration over each interface 
of the domain decomposition method, and a global sum of these local contributions over the structure 

4.4.2 Unstable phase (arc-length procedure) 

The convergence of the algorithm is evaluated by computing Criterion i/^^'^^ after each prediction stage 
of the Newton scheme (the residual of the control equation, which has no physical meaning, is not 
accounted for). 

4.4.3 Typical values 

Our experiments of delamination analysis within the LaTIn framework have shown that a stopping 
criterion i/'*'^'' set to i'^*^'' = 1 x 10~^ ("d" stands here for "desired") is usually sufficient to ensure a 
global convergence of the iterative process (at least, crack fronts are correctly localized, which means 
that the large wavelength piece of information is correctly captured). In our simulations, and in order 
to force an accurate convergence of the local quantities (equilibrium of the interface forces and verifi- 
cation of the cohesive law in the process zones), ly'^'^^ is set to 1 x lO^'^. 



5 AN AUTOMATIC PROCEDURE TO CONTROL THE LOAD 
INCREMENTS 

In this section, we combine the ideas detailed in Section [3] to estimate the time discretization error, 
and the developments of the last section, dedicated to the evaluation of the convergence of the iterative 
parallel resolutions to derive a new time discretization error indicator, suited (but not restricted) to 
the mixed domain decomposition strategy. Based on this new indicator, an automatic procedure to 
control the load increments is derived. 



5.1 Time discretization error criterion in a domain decomposition frame- 
work 

5.1.1 Definition 

A sufficiently converged solution of the reference problem being reached at Current time tn+i, by 
making use of the LaTIn-based resolution strategy, a continuous solution is constructed over [t„ t„+i], 
as described in Section [sj A new time discretization error criterion j/t^e.dd (•"gj^" gtands for "domain 
decomposition") is computed with respect to the interpolated solution: 

time,dd interp,dd /oo\ 

ly,. = max ly.T [62} 

where we recall that Ng + 1 is the number of intermediate times (ii)ig|o, a^,] such that <ti < t„+i 
at which intermediate solutions are constructed, and Criterion i/^^-t^^rp-dd jg evaluated. 
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^interp,dd computed by the same formulas defining v^^ 
variables is known "continuously" over Time interval [0 ti] (from the interpolation): 



except that the history of the interface 



(^1 



q p 



interp.dd^ ^ ^ ^ ^ ^ ^^interp,dd^ wllCrG ^j^interp^dd^ F^F 

j=i j=i 



(33) 



and, in Equation (29), the stiffness operator of the cohesive interfaces ^ ee' {E,E')£li, ueY replaced 



by the reconstructed operator K_ iiW i 



])■ 



Note that Time discretization error criteria ^^^^ ,,time 



-E)\te[0 ii 

^ume,dd ^.^^ ^time g^j,g slightly different. In Section [sj 
we assumed that the solutions obtained at Times t„ and i„+i satisfied the global static admissibility 
(Equation ([2])). This assumption cannot be made anymore if the solver used is the mixed domain 
decomposition strategy (unless the convergence criterion threshold is set to a very low value, which 
would be ineffective, from a numerical point of view). Indeed, the equilibrium is only satisfied in terms 
of substructures and macroscopic interfaces variables. In addition, the kinematic admissibility is not 
fully satisfied, for each ply has been decomposed into substructures which are imperfectly bonded 
during the resolution process. Hence, the new time discretization error criterion i,time,dd jneaguj-gs 
not only the cohesive law residual, but also an interface microscopic equilibrium residual (which is 
small) and a jump of displacement through perfect interfaces, both due to a partial convergence of the 
iterative solver. 




Figure 10: Grey curves : evolution of jy^^^t'^^P^d-d function oft G [tn tn+i] for different values of At. 
Black curve: Evolution of i/*^"^^'<^d ^[^^i respect to At (maximum values of the gray curves) 



5.1.2 Properties 



Figure ( 10 ) shows the evolution of i/'^'^t^^P'dd a time interval [t„ tn+i] for a given computation 

time tn of the unstable dclamination simulation represented Figure ( |12[ ) (which will be detailed later 
on). The different gray curves correspond to different values of the time increment At (value of the 
prescribed arc-length in this case) . Note that the the value of j/«"*erp,dd Computation time t„ is the 
value lyf^^ of v^*^'^ which has been prescribed to ensure a sufficient convergence of the LaTIn iterative 
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process. From this set of parametric studies, the values of be plotted with respect to At 

(maximum values of i/^^-t^^^P^d-d, ^y^^ j^^^ tn+i], black points on the figure). The resulting function (black 
interpolated curve) is monotonic. 

One can also remark that even when a large time step is prescribed, the curve j^^^-tf^^p.dd ^ 
function of i S [i„ tn+i] is smooth. Thus, a relatively small number of evaluation of this residual over 
the time step is sufficient to obtain an accurate value of the time discretization error criterion j/*™^.'*''. 
In addition, as the computation of Criterion i/'^^t'^^P^dd cheap, even a large number of intermediate 
time steps would not alter the numerical efficiency of the strategy. In practice, we choose iVg = 10. 

5.2 Adaptive load increment procedure 

Our aim is to solve the delamination problem at Computation time t„+i under the constraint of given 
value i/'J-'^'^''^'^ of the time discretization error indicator, the current time increment At = tn+i ~ tn 
{i.e.: the prescribed arc-length or the load increment) being set as a new unknown. A quasi-Newton 
technique is used to solve the nonlinear constraint equation: 

^Ume,dd^^^^ _ ^tjme,dd ^ q ^3^^ 

Basically, each iteration of this scheme is decomposed in three steps: 



• a linear step, where a value of the time increment is predicted (see formulas ( 35|36 1 in next 
paragraph). 

• a correction stage where the full delamination problem is solved, at the current time step tn+i, 
until convergence of the underlying nonlinear solver. The time increment At is here fixed to its 
predicted value. 



the computation of a convergence indicator (norm of the residual of Equation (34)) 



The linear prediction stage at Iteration fc + 1 of Computation time consists in solving the 



linearized relation linking jyt^e.dd ^j^g time increment At (see Figure (11 1). This prediction is done 
by the following formula: 

time,dd _ time.dd^^^ 
jytime.dd'^ j^time.dd^ 

Previous formula does not warranty the prediction of a positive arc-length (the function to linearize is 



concave). If a negative time increment is predicted, equation (35) is replaced by the following relation, 
which has empirically shown good convergence properties: 



/ time.dd 

At^+^ ^ J ^ ^ At'' (36) 



6 VALIDATION OF THE STRATEGY 

6.1 First numerical studies of the time step control procedure on the stable 
DCB case 

Let us fully detail the results obtained on the stable DCB case presented in Section [2] Figure [2j 
This problem is globally stable and solved using, at each computation time, the LaTIn-based domain 
decomposition strategy. In the four simulations presented Figure [3j the time increments have been 
obtained by prescribing increasing values ofiyl'""', respectively 1.0 x 10"^, 5.0 x 10"^ 2.0 x 10"! and 
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3.5 X 10~^. The resulting average time increment increases, the total number of computation times N 
being respectively equal to 69, 21, 9 and 5. 

As explained in Section |3j the damage state in the cohesive interfaces tends to the one obtained 
for very small load increments (case 1) when the value of i/^™-'^-'^'^ decreases. More precisely, the 
delaminated area of the cohesive interfaces [i.e.: the dissipated energy) converges in a monotonic 
manner with decreasing values of threshold of the time discretization indicator. When this threshold 
is set to a value smaller than 2.0 x 10"""^, the error made in terms of dissipated energy is not significant. 

Though, this test case is too specific (stable, only one crack front) to give a reliable threshold value 
^time,dd .^j^j^^j^ should bc applied in the general case in order to insure a sufficient convergence of the 
solution with respect to time. 



6.2 Unstable holed-plate delamination problem 




Figure 12: Definition of the holed plate problem (317 000 d.o.f.) 



We consider a eight-plies holed-plate structure, under traction (prescribed displacements). The 
plies are orthotropic (stiffness ratio: 1/20) and the stacking sequence is [0 ±45 90] s, which leads to the 
initiation of delamination due to edge effects. The initial stiffness properties of the cohesive interfaces 
are obtained by the same homogenization in the "thickness" of the interfaces (one tenth of the thickness 
of the plies) which has been described for the DCB problem in Section [2] Due to the material and 
structural symmetries, only the top half of the structure is simulated. The unstable quasi-static time 
analysis is performed by making use of the arc-length procedure described in Section [4] The global 
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response curve (plotted in Figure (13), Case 3) of this case shows two main zones of instabiUty. The 
first one corresponds to a an unstable propagation of the delamination in the [—45 + 45] interface 
while the second one is a crack propagation in the [0 — 45] interface, both mainly in shear mode. 



6.2.1 Prescribed time step (Cases numbered 1, 2 and 3 in the whole analysis of the 
results) 

The first set of simulations is performed by successively prescribing three different fixed arc-lengths. 
The arc-length which has been arbitrarily chosen in Case 1 is divided by three in Case 2, and by nine 



in Case 3. Instabilities appear in the global response of the structure (Figure (13)). Figure (14) shows 
the damage maps in the [0 — 45] [—45 + 45] and [—45 90] cohesive interfaces in a monotonic phase 
of the global behavior (limit point after which the delamination front evolves in an unstable manner 
in the [0 — 45] interface, which corresponds to the circled point in all graphs of Figure (13)). This 
particular point of interest has been reached respectively in 8, 63 and 237 time increments. Note that 
we do not aim at discussing the validity of the solutions reached but at ensuring that the incremental 
strategy follows the equilibrium path of the converged solution with respect to the time. 
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Figure 13; Global reaction force versus prescribed displacement in the holed plate problem under 
traction, three different predefined arc-lengths being applied (Cases 1,2 and 3) 

No significant difference can be observed in the damage maps and global response curves corre- 
sponding to the two finest analysis in time, which means that the solutions are sufficiently converged 
with respect to time in Cases 2 and 3. In Case 1, the time increments are too coarse, which results 
in the incremental resolution procedure to follow a different equilibrium path (see the damage maps 



in Figure (14)). This phenomenon can impair the global response of the structure, as it can be seen 
on Figure (13). The instability phases framed on the converged solutions (Cases 2 and 3) are wrongly 
predicted in Case 1. These differences appear even more clearly on the dissipated energy versus pre- 
scribed displacement curves plotted on Figure ( [T5| (the curves labeled "reference" and "coarse grid" 
refer respectively to Cases 3 and 1), corresponding to the first global instability and to the following 
stable phases of the time analysis. 



Figure (|16|) presents the values of (j^|*(™'^'°''^)neli, JV] a function of the successive computation 



times in Cases 1 and 2, from the beginning of the analysis to the starting point of the second global 
instability. One can see that in Case 2, in which the time increments are sufficiently small to let the 
iterative algorithm follow the correct equilibrium path, the values of the discretization error indicator 
from 1 X 10 to 1 X 10 ^. Conversely, we show in the next set of studies that setting the 
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Figure 14: Damage state in the interfaces of the holed plate at the beginning of a global instability 
in the case of a coarse time grid (Case 1), and in a converged case (Case 2). A Fixed arc-length is 
prescribed in both cases. In the first case, the damage in the [—45 + 45] interface is underestimated. 



threshold value jy**™^- of the time control procedure to the maximum of the values iit^™'^^d.d obtained 
in Case 2 permits to obtain a correctly predicted solution. 



6.2.2 Control of the time step (Cases numbered 4 and 5) 

The second set of simulations is performed by making use of the procedure described in Section [5] 
to control the successive prescribed arc-length, j^*™'^-'*'' jg successively set to 1 x 10~^ (Case 4) and 
5 X 10^'^ (Case 5). The damage state in the cohesive interfaces at the beginning of the second global 
instability phase predicted by prescribing j^**™'^''''* = 1 x 10~^ is very closed to the one obtained in 



Case 2 of our first set of simulations (see Figure (17)). The total number of time increments drops 
from 63 to 40. 

As explained previously, the dissipated energy versus prescribed displacement curves (Figure ( [T5| ) 
obtained in the reference case 3 (very small prescribed arc-length) and in Test case 1 (coarse prescribed 
arc- length) are very different (incorrect equilibrium path in the second case). When using the time 
control strategy, i/^'^'^''^'^ being successively set to 1 x 10^^ and 5 x 10^'^, the correct equilibrium path 
is followed. In addition, the dissipated energy versus prescribed displacement curves gets closer to the 
one obtained in the reference simulation when the value of j/*™*^"'*'^ decreases. 



6.3 Discussion 

The threshold value found numerically here can be compared to the one prescribed to ensure the 
convergence of the iterative resolution strategy at each computation time, i'^^'^^ ■ As explained in 
Section |4j the value v]^'^'^ which ensures a sufficient convergence of the LaTIn solver can be obtained 
empirically by performing time independent benchmark tests (for instance the first time step of a 
delamination analysis). The time control strategy developed in this paper consists in monitoring the 
residual of the reference problem equations continuously during the time analysis, the measure used 
at any time being a time independent version of Hence, it is not surprising to find out in the 

numerical examples that the higher value jy^J-^^''^'^ permitting to follow the correct equilibrium path is 
the value of i/**'^'' which permits to obtain the convergence of the global informations (position of the 
crack fronts) at each computation time. Hence, applying the time control procedure only requires the 
prior knowledge of indicator v^J^"^^ ■ 
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prescribed displacement (mm) 



Figure 15: Dissipation versus loading curves for different resolution strategies: explicit fine and coarse 
time steps (Cases 1 and 3, fixed arc-length) or automatically controlled time increments (Cases 4 and 
5, fixed time discretization error) 

7 CONCLUSION 

In this paper, we presented a strategy to adapt automatically the time increment in quasi-static de- 
lamination problems to the very sharp non-linearities which are involved. This strategy is based a 
continuous monitoring of the residual of the reference problem equations with respect to time. This 
has been achieved by calculating the error in the constitutive law on admissible solutions interpolated 
over each time steps, which enables to define a time discretization error criterion evaluating the rel- 
evancy of the nonlinear computations performed at each time increment. Based on this indicator, a 
simple procedure to control the time step has been derived. The main parameter of this technique is 
easy to obtain as it only requires to perform time-independent benchmark tests prior to the delami- 
nation simulations. The validity of this procedure has been demonstrated on delamination problems 
undergoing global instabilities. 

Our current interest being to perform buckling-driven delamination analysis, the validity of this 
strategy shall be verified, in the future, on computations involving geometrical nonlinearities. 
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